clear all; clc;
syms x h;

M = inv([1 0 0 0; 0 1 0 0; 1 h h^2 h^3; 0 1 2*h 3*h^2]);

for j = 1:4,
    a = M(:,j);
    phi(j) = sum([1 x x^2 x^3]*a);    
    d2phi(j) = diff(phi(j),x,2);
end

disp('Using basis functions that have u as expansion coefficient:')
disp('A(i,i):')
pretty(int(d2phi(1)*d2phi(1),x,0,h)+int(d2phi(3)*d2phi(3),x,0,h))
disp('A(i,i-1):')
pretty(int(d2phi(1)*d2phi(3),x,0,h))

disp('Using basis functions that have du as expansion coefficient:')
disp('A(i,i):')
pretty(int(d2phi(2)*d2phi(2),x,0,h)+int(d2phi(4)*d2phi(4),x,0,h))
disp('A(i,i-1):')
pretty(int(d2phi(2)*d2phi(4),x,0,h))

disp('Using a mix of basis functions that have u or du as expansion coefficient:')
disp('A(i,i) (two possible combinations ensure continuity, both with the same value):')
pretty(int(d2phi(3)*d2phi(4),x,0,h)+int(d2phi(1)*d2phi(2),x,0,h))
disp('A(i,i-1) (two possible combinations ensure continuity, both with the same value):')
pretty(int(d2phi(2)*d2phi(3),x,0,h))

disp('Since only load is a point load at the end, f(x) equals zero everyehere expect for direc delta function at the end, picking out the value of phi(j) at that node')